## Replicate Reanalysis of King, et al. (2007) design
## Paper Figure 1
## 3 May 2012

## Covars for RUR, URB
covars.rur.match <- c("pobtot",  "gradoes", "pro_ocvp","percfem", "perc0_4",
                      "perc18",  "sindss",  "condss",  "res95",   "condisc", "analf15",
                      "casada",  "catolic", "ocupada", "ocusecp", "ocusecs", "ocusect",
                      "perc_2sm","hogjeff", "p5hli",   "margw",   "camas_no","consu",  
                      "medic",   "enferm",  "inc_spss","superfi", "pobno1k", "ppardes",
                      "ptecdes", "ppisdes", "pelectr", "prefri",  "paguent", "pdrenaj",
                      "pdreagu", "percop",  "infra120")

covars.urb.match <- c("pobtot",  "gradoes", "pro_ocvp","percfem",
                      "perc0_4", "perc18",  "perc60",  "perc65",  "ppsderss","ppcderss",
                      "ppimss",  "res95",   "condisc", "alf15",   "casada",  "catolic",
                      "ocupada", "ocusecs", "ocusect", "perc_2sm","phogjeff","margw",
                      "camcens",  "camnoce",  "quiro",    "consu",    "medic",  "enferm",
                      "altitud", "ppardes", "ptecdes",  "ppisdes",  "pelectr",  "prefri",
                      "paguent",  "pdrenaj", "pdreagu",  "percop")

load("repRurdata.mmgosp.comb.RData")
load("repRurvc1.RData")
load("repUrbmmgosj1.RData")
load("repUrbvc1.RData")

library(blockTools)
rurblock <- block(data = rurdata.mmgosp.comb, vcov.data = rurvc1, groups =
  "cve_edo", id.vars = "clues", block.vars =
  covars.rur.match)
urbblock <- block(data = urbmmgosj1, vcov.data = urbvc1, groups =
  "claveedo", id.vars = "clues", block.vars =
  covars.urb.match)

## Cut unvalidated pairs
source("unvalCut.R")
rurblock <- unval.rurcut(rurblock)
urbblock <- unval.urbcut(urbblock)

rurblocknaive <- block(data = rurdata.mmgosp.comb, vcov.data = rurvc1, groups =
  "cve_edo", id.vars = "clues", block.vars =
  covars.rur.match, algorithm="naiveGreedy")
rurblocknaive <- unval.rurcut(rurblocknaive)

for(i in 1:length(rurblock$blocks)){
  o <- paste(sum(rurblock$blocks[[i]][,3] <=
    rurblocknaive$blocks[[i]][,3]), " of ",
             nrow(rurblock$blocks[[i]]), " better/same under optGreedy, state number ", names(rurblock$blocks)[i], sep="")
  print(o)
}

urbblocknaive <- block(data = urbmmgosj1, vcov.data = urbvc1, groups =
  "claveedo", id.vars = "clues", block.vars =
  covars.urb.match, algorith="naiveGreedy")
urbblocknaive <- unval.urbcut(urbblocknaive)
rm(unval.rurcut, unval.urbcut)

for(i in 1:length(urbblock$blocks)){
  o <- paste(sum(urbblock$blocks[[i]][,3] <=
    urbblocknaive$blocks[[i]][,3]), " of ",
             nrow(urbblock$blocks[[i]]), " better/same under optGreedy, state number ", names(urbblock$blocks)[i], sep="")
  print(o)
}
print(paste(1+31+12+2+2+6+1+1+8+1, " better/same under optGreedy.", sep="")) ##[1] 65
print(paste(1+35+12+3+2+2+6+1+1+9+1+1, " total.", sep="")) ##[1] 74

storage <- expand.grid(meas = c("mean", "median"), rural = c("rural", "urban"), alg = c("optgreedy", "naive"), set = c("all", "half"), dist = NA)

mat <- rurblock$blocks[[1]]
for(i in 2:6){
  mat <- rbind(mat, rurblock$blocks[[i]])
}
storage$dist[1] <- mean(mat[,3]) ## r 5.05
storage$dist[2] <- median(mat[,3]) ## r 4.229949
mat <- urbblock$blocks[[1]]
for(i in 2:6){
  mat <- rbind(mat, urbblock$blocks[[i]])
}
storage$dist[3] <- mean(mat[,3]) ## u 5.38385
storage$dist[4] <- median(mat[,3]) ##  u 4.657221

matn <- rurblocknaive$blocks[[1]]
for(i in 2:6){
  matn <- rbind(matn, rurblocknaive$blocks[[i]])
}   
storage$dist[5] <- mean(matn[,3])  ## r 5.08
storage$dist[6] <- median(matn[,3]) ## r 4.832526
matn <- urbblocknaive$blocks[[1]]
for(i in 2:6){
  matn <- rbind(matn, urbblocknaive$blocks[[i]])
}   
storage$dist[7] <- mean(matn[,3])  ##  u 5.955974
storage$dist[8] <- median(matn[,3]) ##  u 5.439292

mat <- rurblock$blocks[[1]][1:ceiling(nrow(rurblock$blocks[[1]])/2),]
for(i in 2:6){
  mat <- rbind(mat,
               rurblock$blocks[[i]][1:ceiling(nrow(rurblock$blocks[[i]])/2),]) 
}
storage$dist[9] <- mean(mat[,3])  ##[1] r 3.663548
storage$dist[10] <- median(mat[,3])  ##[1] r 3.315963
mat <- urbblock$blocks[[1]][1:ceiling(nrow(urbblock$blocks[[1]])/2),]
for(i in 2:6){
  mat <- rbind(mat,
               urbblock$blocks[[i]][1:ceiling(nrow(urbblock$blocks[[i]])/2),]) 
}
storage$dist[11] <- mean(mat[,3])  ##[1] u 5.170947
storage$dist[12] <- median(mat[,3])  ##[1] u 4.068457

matn <- rurblocknaive$blocks[[1]][1:ceiling(nrow(rurblocknaive$blocks[[1]])/2),]
for(i in 2:6){
  matn <- rbind(matn, rurblocknaive$blocks[[i]][1:ceiling(nrow(rurblocknaive$blocks[[i]])/2),])
}   
storage$dist[13] <- mean(matn[,3])  ## r 3.958523
storage$dist[14] <- median(matn[,3]) ## r 3.803211
matn <- urbblocknaive$blocks[[1]][1:ceiling(nrow(urbblocknaive$blocks[[1]])/2),]
for(i in 2:6){
  matn <- rbind(matn, urbblocknaive$blocks[[i]][1:ceiling(nrow(urbblocknaive$blocks[[i]])/2),])
}   
storage$dist[15] <- mean(matn[,3])  ## u 5.833643
storage$dist[16] <- median(matn[,3]) ## u 4.668452
